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Abstract 

Time-continuous dynamical systems denned on graphs are often used to model complex 
systems with many interacting components in a non-spatial context. In the reverse sense 
attaching meaningful dynamics to given 'interaction diagrams' is a central bottleneck prob- 
lem in many application areas, especially in cell biology where various such diagrams with 
different conventions describing molecular regulation are presently in use. In most situa- 
tions these diagrams can only be interpreted by the use of both discrete and continuous 
variables during the modelling process, corresponding to both deterministic and stochastic 
hybrid dynamics. The conventions in genetics are well-known, and therefore we use this 
field for illustration purposes. In [25] and [26] the authors showed that with the help of a 
multi-scale analysis stochastic systems with both continuous variables and finite state spaces 
can be approximated by dynamical systems whose leading order time evolution is given by 
a combination of ordinary differential equations (ODEs) and Markov chains. The leading 
order term in these dynamical systems is called average dynamics and turns out to be an 
adequate concept to analyse a class of simplified hybrid systems. Once the dynamics is de- 
fined the mutual interaction of both ODEs and Markov chains can be analysed through the 
(reverse) introduction of the so called Interaction Graph, a concept originally invented for 
time-continuous dynamical systems, see [5]. Here we transfer this graph concept to the aver- 
age dynamics, which itself is introduced as an heuristic tool to construct models of reaction 
or contact networks. The graphical concepts introduced form the basis for any subsequent 
study of the qualitative properties of hybrid models in terms of connectivity and (feedback) 
loop formation. 

1 Introduction 

The modelling of systems with many interacting components usually involves employing differ- 
ent mathematical methods and approaches in order to obtain effective descriptions. A typical 
such situation is the study of cellular systems involving different interactions of (bio-) molecular 
species establishing some cellular function. In this class of problems many processes occur simul- 
taneously, and they are so interconnected that the concept of a network turns out to be a natural 
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framework to work with. It has become standard in cell biology to draw diagrams to illustrate 
the molecular systems, where only part of the components and arrows can be straightforwardly 
interpreted as classical chemical reactions. In order to introduce a well-defined network struc- 
ture it is necessary to define a set of interpretation rules relating the system components to the 
network elements. The latter are given by the vertices and edges of the network (compare this 
with the introduction in [E]). There is however no unique such mapping even if the network 
components are well defined. This has been for example illustrated in [5] where different graphs 
are shown that are associated to mass-action reaction networks. But nevertheless there is a quite 
generic such network structure given by the so-called interaction graph X (denoted as Gj in [5] 
). For the setting in this paper the classical interaction graph X needs to be extended, as we 
are now looking at situations where some of the players or particles are not abundant enough to 
be described by a continuous variable. The vertex set (nodes) becomes larger, there is a node 
set containing components describing a single continuum variable (like a chemical concentration, 
which in this framework are the macro-states modeled with deterministic dynamics), and nodes 
associated to every discrete state of the Markov chain describing microscopic interactions. The 
approach taken in this article is purely graphical, i.e. we define vertex and edge sets of graphs 
associated to hybrid dynamics and their scalings. We emphasise that alternatively one could 
introduce the respective adjacency and incidence matrices and present all graph concepts intro- 
duced on an pure algebraic level, [3j[10]. Especially purely algebraic conditions for the uniqueness 
of an invariant measure of a Markov chain will be much more familiar to most readers, see section 
|2.3| where we use graph connectivity as a concept instead. Knowing about this mathematical 
duality we therefore restrict ourselves here to the graph aspect only in order to (a) emphasise 
the connection to (dynamical) network theory, and (b) showing the connection to interpretation 
of (biological) graphical diagrams in a complete and consistent way. 

In this paper we always assume spatial homogeneity. We are interested in a special class of 
systems which are characterised by the presence of two different kinds of variables. One variable 
subset admits a description through the notion of continuity (or concentration) based on the 
assumption of mass action kinetics, which itself is based implicitly by a so-called continuum 
limit. Therefore ordinary differential equations (ODE) are used to describe their dynamics. The 
other variables can only take finitely many possible states and are inherently stochastic. Hence 
the overall system is typically hybrid. Hybrid systems have been used to describe many different 
situations, see for instance [4]. Very close to our framework are [13], and the late work by 
Gillespie [9] . Note that in our framework the dynamics of both the discrete and the continuous 
state spaces have in principle a similar mathematical description. In fact, any set of reactions 
has to be described via stochastic dynamics, typically a master equation (ME) giving the time 
evolution of the probability of the system while being in a certain state at time t. As described in 
[23 ESI EZl HH] microscopically the state of a reaction or 'rule-based' system is an element in x E, 
where is a lattice, S is its scale - possibly a vector as particle properties may scale differently 
for each species - and E is a finite set, the discrete state space. The 'continuous' state space 
accommodates the states of the reactions or rules that can be described by mass-action kinetics 
in a large number limit, whereas E contains the states that are modelled intrinsically discrete. 
Typical examples are conformational changes and binding/unbinding events of macro- molecules, 
but some switch in an engineering example is leading towards the same kind of description. 
In a larger system there are in general many molecular subsystems whose states are finite and 
which participate in cellular regulation. In a cell biology context we refer to such subsystems 
as molecular machines. The consequence of this is that E collects the possible states of all the 
molecular machines and hence E is naturally given by the Cartesian product of the discrete 
state spaces of each such subsystem. In graph terms this generally means that the Markov 
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chain transition graph has as many disconnected subgraphs (internally connected components) 
as there are different such molecular machines. The influence of the algebraic structure of E on 
the dynamics has been partially studied in [26]. In [25] the continuum limit 

h 5 R N as S 

was introduced. In such a regime the occupation numbers in become concentrations x in M. N , 
but the states in E retain their discrete structure. The ME becomes a vectorial Fokker-Planck 
equation of the form: 



= J2 CIA*] ■ M*,x) + - e Y / *w(x) • M*,x) Va e S. (1) 
Here p a (£,x) is a probability density, 

/ dx V/) a (t,x) = 1, 

and e is a time scale describing the ratio of speed of changes in the deterministic part - 
compared to the switching of states in the discrete part of the state space described by the 
Markov chain. The differential operator £[x] is matrix- valued and describes the reactions of mass- 
action type, whereas /C(x) is a Markov chain (MC) generator which describes the interactions 
in E. The parameter e is a scaling factor involved in the continuum limit and used to derive 
the differential operator. Such a parameter is usually infinitesimal and hence ([!]) becomes a 
perturbation problem as e — >• 0. It is not difficult to see that e — >• corresponds to having a 
processes in E which is much faster than a process driven by C\x\. Therefore (111) can be studied by 
an expansion in e through the so called adiabatic theory (see [25] [26] El] ) • The solution for e = 
corresponds to a system staying in the MC equilibrium state, which is an invariant probability 
measure. For e > sufficiently small the system evolves according to a vector field which is an 
average of all possible dynamics indexed by E, and taken against the invariant measure. The 
leading order term of the e-expansion of the solution of is a deterministic equation for the 
probability density (Liouville Equation). This corresponds to a deterministic vector field in the 
concentration space R^. We termed such a vector field the average dynamics. This name for the 
leading part in the expansion derived by the adiabatic theory has been suggested by its derivation: 
the average dynamics is constructed out of the invariant measure /i(x) of the MC together with 
the deterministic dynamics given for each state of the MC. Let /i(x) be the invariant measure 
of the MC and let f( a \x.) be the deterministic vector field for each a G E. Then the average 
dynamics is given by averaging the vector field against the measure /i(x) yielding 

*[x] = £Mx)/ (<T) (x). (2) 

The whole scaling process can be visualized through the diagram in Figure [T] In this paper we 
want to consider the averaging principle only as a heuristic tool to construct meaningful models 
for applications. In this respect we present a method to derive deterministic dynamics from hy- 
brid stochastic dynamics. We term this modeling procedure the averaging principle. We further 
assume that the averaging principle applies to the specific system (which means can be modeled 
correctly) , but also include an example where some unstructured particles do not form a density 
at the end of the paper. The deterministic vector field is naturally associated to an interaction 
graph and similarly another interaction graph is associated to the MC (denoted by Xmc)- 
These two graphs are combined into a larger graph Xc describing how stochastic and determin- 
istic variables interact. The average dynamics describes the time evolution of the deterministic 
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Figure 1: A system is partitioned into two "subsystems". One subsystem is driven by determinis- 
tic interactions, which are represented by a classical interaction graph. The other subsystem has 
a stochastic evolution dictated by a MC that has a natural interaction graph interpretation. The 
two subsystem are coupled by a stochastic dynamics which is described only at the level of the 
Master Equation. The continuum limit and the adiabatic theory lead to the average dynamics, 
the object is studied in this paper. 
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part averaged against the invariant probability measure. In fact the invariant measure assigns 
probabilistic "weights" to the finite states of the MC. To each such state corresponds a specific 
deterministic vector field which will give a contribution to the average vector field according to 
the probabilistic weight of the MC state. The new vector field produced by the average proce- 
dure is then associated to a new (in a sense collapsed) interaction graph X. The new collapsed 
graph resonates with the notion of condensation of a graph [31 . In X the subgraph Xmc is no 
longer present but new edges will appear. These new links connect vertices which were part of 
the Markov chain graph Xmc before applying the averaging procedure (see Figure [l]) . Similar 
models are already known in the literature (see [21]), but the average dynamics presented here is 
derived from basic principles that correspond (asymptotically) to the deterministic limit of the 
originally stochastic system. 

The aim of this paper is to introduce and explain all the potentiality of the averaging principle 
and the associated averaged interaction graph X. We first describe the notion of the interaction 
graph for deterministic dynamics and for Markov chains separately, and then extend the concepts 
to hybrid and averaged systems. Furthermore we investigate the consequences of the structure 
of E from the modeling point of view. The space E is in general a Cartesian product of many 
discrete state spaces. We show that the structure of the infinitesimal generator JC determines 
whether the discrete state (molecular) machines behave as independent subsystems. Finally 
the averaging principle is illustrated by several examples taken from genetics. In particular we 
shall show that the dynamics of so-called feed-forward loops as analysed in [20] can be derived 
more rigorously (in the sense that the underlying scalings to derive the macroscopic equations 
are at least known) by the average dynamics. As a final example we show that the average 
dynamics can be used to re- formulate a model for genetic oscillators as introduced in [30]. This 
example shows that oscillations are produced by the interaction of two genes. We show that 
the corresponding MC can effectively often be decomposed/reduced into two independent MC's 
describing the dynamics of each single gene. The evidence of the rhythms for this model is shown 
numerically. The modeling techniques described in this paper give also rise to model existing 
synthetic genomes from a dynamical perspective, with the hope to achieve that in future in a 
semi-automatic way once a sequence is known, see [15] [16] . 

2 Interaction graphs 

2.1 The interaction graph for deterministic dynamics 

Let us consider a dynamical system defined by N first-order ordinary differential equations of 
the form 

Xi{t) = . . . , x N (t), a) + 9i (t, 0), i = 1, N. (3) 

Here x = (xi, ..,xn) G is the continuous state of the system, a, /3 are (sets of) parameters, 
and / = (/i,..,/jv) and g = (gi,..,#jv) are both sufficiently smooth vector fields on R N . We 
associate a directed graph to which we call the interaction graph X of system If g is not 
identically zero then (which implies ([3| is non-autonomous) an additional state variable xq is 
introduced. This state xo (not affected by other states) is then called the environment. 

Definition 2.1 (The interaction graph Xp for deterministic dynamics). The interaction graph 
associated to tew is a directed graph denoted by Xd = (V, E) where 
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(i) the vertex set V(Td) is equal to the collection {xi, ...,xn} of state variables in case g = 
in and {xojXij ...jXn} otherwise, 

(ii) an edge = (xi, xj), with z, j = 1, . . . , N is an element of E(Xd) iff 

dfiU) 

— -p- is not identically zero. 

OXi 

An edge e^j = (xq, Xj), with j — 1, . . . , N is an element of E iff gi is not identically zero. 
(Hi) The edge is directed from xi to Xj. This is equivalent to say (x^Xj) is an ordered pair. 

Remark 2.1. The vertex xq, the environment, describes any external actions/influences on the 
system. Note that in applications external actions may be used to model a dynamics which is 
much slower than the one of the system and consequently such external variables can be considered 
"frozen". Obviously the edges of Id can be weighted with the sign of telling whether an 

interaction is positive or negative. If such sign patterns are identical for all states x of the 
respective state space, then the system is called quasi-monotone and more tools are available to 
characterise the qualitative properties of the system fWj/. The definition ofXjj can be generalised 
by including time dependent vector fields, that would imply a time dependence of the link structure 
in Xjj as well. But even for general autonomous ODE generating a flow the state is time- 
dependent, so Xjj does possibly change its links/edges and weights. This of course does not hold 
for initial conditions identical to an equilibrium. 

Remark 2.2. Note that using the vector field we can readily construct the adjacency matrix A 
ofXrj by defining: 



otherwise 



2.2 The interaction graph for Markov chains 

A linear Markov chain (MC) has a natural structure of a directed graph which we denote by 
Xmc- 

Definition 2.2 (The interaction graph Xmc fc> r Markov chain dynamics). (i) V(Xmc) — {o~i-, • • , cfm}, 
where E = {<T l5 . . . , o~m} is the collection of (finite) MC states, 

(ii) for each couple of vertices G £ x E an edge is an element of E (Xmc) iff the 

probability of the transition i —> j is different from zero. The edge is directed from i to 
j- 

Remark 2.3. Note that in general transition probabilities are functions of x and therefore in 
general time- dependent. This implies again the edge set E of Xmc may change in time. 

In many situations it is more useful to describe a MC through its infinitesimal generators rather 
than by the transition probabilities. The infinitesimal generator of a MC is determined by the 
rates at which transitions occur. The rates themselves correspond to probabilities restricted 
to an infinitesimal time interval (see [29]). Let E = (<ji, cfm) be the states of the MC un- 
der consideration. At time t the state of the MC is determined according to the probability 
distribution 
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P{t) = (PiW,.--,PmW), 

where Pi(t) is the probability that at time t the system is in state <r^. The respective normalisation 
condition is given by 

M 



!>(*) = L 



Given two states z, j, the transition i —> j occurring in a time interval St is determined by the 
probability 



P 



/Qj5£ for i ^ j, 



where /Qj is the f/-th element of the infinitesimal generator. Letting St —> the infinitesimal 
form of the dynamics can be written in matrix form: 

M 

Pi(*) = X)ft'W /C ^ i = l,-..,M. (4) 

Remark 2.4. In general infinitesimal generators depend on various parameters which might 
have their own dynamics and time dependence. We are interested in autonomous systems where 
the dynamics associated to the Markov chain transitions are all faster than the rest of the ( in this 
paper mostly deterministic) dynamics which then may affect JC on a longer time scale. Therefore 
we assume JC not to be explicitly dependent on time, only implicitly via a possible dependence on x. 

We are now able to state an alternative definition of Imc in terms of the infinitesimal generator 
JC: 

Definition 2.3 (The interaction graph Imc for Markov chain dynamics based on JC). Given 
the MC interaction graph is a directed graph denoted by Imc = (£,£'), with V(1mc) = 
{<ti, . . . , o~ M } as before, and 

(i) an edge (arrow) e*ij is an element of E iff the couple of vertices (i,j) (with i ^ j) is such 
that JCij is not identically zero. Moreover for every i G S ; en is not in the set E. 

(ii) the edge e^ is directed from i to j and is also denoted by i — » j . 

Remark 2.5. From the normalisation condition one can show (see [29]) that the infinitesimal 
generator satisfies the following condition: 

M 

^ JC^ = for every i = 1, . . . , M. (5) 

j=i 

Remark 2.6. Note that using JC for Imc we can readily construct the adjacency matrix A by 
defining: 

A _ j * if JCij ^ 



{i 



otherwise 
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Condition ^ implies that JC has a non-trivial left and right kernel. This is one of the main 
properties characterising a MC. 

Remark 2.7. It holds that \i = (/ii, ..,/xm) is a stationary measure if 

M 

^] ^jJCji = for every i = 1, . . . , M, (6) 

and 

M 

If /ij > for all j , then \i is said to be strictly positive. 

For the spectrum of JC we recall the following result: 
Lemma 2.1. The spectrum of JC is contained in the left half of the complex plane. 

This result can be proved by using Gershgorin's theorem (see [22]). 

2.2.1 Convex combinations of stationary measures 

Equation ^ defines the left kernel of the matrix JC (or the kernel of its transpose JC T ). We 
denote the left kernel of JC by ker(/C) and assume that dim(ker(/C)) = L > 1, namely ker(/C) is 
generated by L stationary measures {/i/}^. Under this condition we define 

Definition 2.4 (Convex combinations of stationary measures). Let be the stationary 

measures generating ker(/C). £e£ {0/}^ =1 be L positive real numbers with ^2f = i@i — 1. W^e caZZ 
t/ie following sum 

tte convex combination of stationary measures associated to ker(/C). 
The convex combinations have the following properties: 

1. Any convex combination is a probability measure. 

2. Any convex combination is a stationary measure for the MC. 

Stationary measures are essential to describe the asymptotic behaviour of the MC. In fact one 
can show 

Theorem 2.1. Given any initial probability measure po = (pi(0), ...,pm(0)) ; there exist {Oi}f =1 
with Y^d=i ®i ~ 1 such that the solution p(t) of with p(0) = po satisfies 

lim p(t) = fi 

with 

z=l 

TTie measure n is called invariant measure. 
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Proof. This theorem rephrases some known results that can be found for example in Stroock 



• The matrix /C generates a transition probability matrix by 

p(t) = fi P(t) where P(t) = exp(/Ct), 
where /io is the initial probability distribution. 

• Invariant measures \i are left-eigenvector of P with eigenvalue 1 and /i's form a convex set. 

• The space £ can partitioned into equivalence classes of states 

C = {i, j e E : P i:j (t) > 0, P^(t) > for some t > 0}. 

• The limit lim^ +00 p(£) converges to an invariant measure supported on the union of the 
equivalence classes C. 

If the number of equivalence classes is L a general invariant measure can be written as 

1 = 1 

□ 

The preceding result shows that the asymptotic behaviour for t — >• +oo of a Markov chain is 
classified by the stationary measures, namely by the solutions of equation |6|. In the next section 
we shall describe the connections between the ergodic properties and the interaction graph in 
Markov chains. 

2.3 Interaction graph X M c and some ergodic properties 

In this paragraph we consider the ergodic properties of one single indecomposable Markov chain 
in terms of connectivity of its associated interaction graph. This indecomposable Markov chain 
will be associated to one single molecular machine being present. In case several such machines 
are present (as discussed later) the Cartesian product of state spaces gives rise to several discon- 
nected components or subgraphs (which themselves of course should be internally connected in 
order to represent a meaningful single molecular machine). Such subgraphs can be discussed in- 
dependently of each other, so without loss of generality we assume in this section there is only one 
connected component present. Let therefore (E,/C) be a Markov chain with (an unstructured) 
state space E and generator /C. Let us define 

Definition 2.5. A directed path i j in Imc is a sequence {ii, ...i n } C E such that 

i -> h -> i 2 ... -> i n j- 

Definition 2.6. A vertex i is said to be connected to a vertex j if i =$> j. 

In the theory of directed graphs (see for example [8 ) there are two possible definitions of 
connectivity. Let E be a non-empty set and consider the directed graph (£,22): 

(£,22) is a weakly connected graph if for every z,j G E at least one of the paths i => j, 
j i is not empty. 
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(E,i£) is a strongly connected graph if for every i,j G E both the paths i =>• j , j =>• i are 
not empty. 

One can observe that a strongly connected Xmc represents a MC with a unique invariant 
measure /i, with Hi > for all i G E. A weakly connected Xmc corresponds to a MC with a 
unique /i, but with possibly some zero components. The reason for this outcome is that invariant 
measures tend to concentrate in certain vertices which are ends of directed paths. This motivates 
the definition of end-vertices: 

Definition 2.7. A vertex i e G E is said to be an end- vertex if 

• there exists i G E such that i => i e , 

• and there is no j G E such that i e => j. 

Remark 2.8. In the context of Markov chain theory i e is an absorbing state (see [SI]). An 
absorbing Markov chain contains at least one absorbing state (see f3D]). 

Now we propose to look at connected Xmc which have at most one end- vertex. Such graphs 
will correspond to Markov chains with unique invariant measure that need not be necessarily 
strictly positive: 

Definition 2.8. The graph Xmc is called ^-connected if 

• Xmc is connected, 

• Xmc contains at most one end-vertex (absorbing state). 

Remark 2.9. Note that if a graph is strongly connected then it is ^-connected. If it is 
connected, then it is necessarily weakly connected. 

Proposition 2.1. 7/(E,/C) is such that Xmc is *- connected, then (E,/C) has a unique invariant 
measure. 

Proof. It suffices to show that a ^-connected E with M states has a generator JC with rank(/C) = 
M — 1. In the proof we shall construct interaction graphs by adding a vertex to a given graph. 
This process must be done without altering the property of being ^-connected, i.e. without 
increasing the number of end- vertices. We proceed by induction on M. For M = 2, the graph is 

1 -> 2, 

with generator 

1^ ( ku fcio \ n 7 7 
JC 2 = ( Q I and kn = -k 12 . 

Then one can readily show that 

rank(/C2) = 1. 

The addition of a vertex leads to 

1 -> 2 -> 3 
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/C 3 = I k 22 k 23 , with kn = -k 12 , k 2 2 = -k 23 . 
\ / 

The rank of JC3 is equal to 2 since rank(/C2) = 1, i.e. k± 2 ^ 0. Consider (Sm-i,^m-i), with 
M — 1 vertices and rank(ZCM-i) = M — 2. The operation of adding a vertex leads to a new Sm 
with M vertices and a generator 

JCm = ^ ^ ^ 5 w ith A: a column vector of dimension M — 1. 

The vector /c gives the rates of the transitions from the vertices in Sm-i to the added M th 
vertex, hence k is different from the zero vector. Up to re-numbering the vertices we can consider 
in Em the vertex M — 1 linked to the vertex M, hence &m-i m/0 and necessarily kM-i m-i 7^ 0; 
therefore det(ZCM-i) 7^ yields rank(/CM) = M — 1. □ 

Note that there are stronger purely algebraic conditions for the uniqueness of an invariant 
measure. The previous 'graphical' result can also be shown to be a simple consequence of 

Proposition 2.2. A is an eigenvalue of JC if and only if it is an eigenvalue of a submatrix 
associated to a strongly connected component of the interaction graph associated to K,. 

For the well-known proof of the proposition see for example [31 . For ^-connected graphs 
the unique end- vertex forms itself a strongly connected component that yields a zero eigenvalue 
which is therefore unique, and so leading to rank(/CM) = M — 1. 

The two interaction graphs Xd and Tmc introduced so far can be interpreted as describing 
two different scales of the processes which now need to be related in a further modelling step. The 
deterministic interaction described by Xd is typically describing macroscopic properties, such as 
the overall 'systems' concentrations of different species of very abundant small molecules. Often 
in modelling these continuous variables forming the nodes of Tp will be used for 'communication' 
between the states of the machinery described by the Markov chain. In many cases, such as 
in stirred chemical tank reactors there is no need to introduce any further discrete states for 
the system at all. The process description can be closed by introducing a set of all possible 
chemical reactions and - by using mass action kinetics - to derive the respective time-continuous 
dynamical system describing the chemical reactor macroscopically. The states of the Markov 
chain are connected to 'system switches' of one or several machines govering the processes. In 
our examples in the current paper such machines will be genes. They do not form concentrations 
of any kind that would need to be modelled in order to close the model description. 



3 The averaging principle, combined and averaged inter- 
action graphs 

We now introduce the averaging principle which has been derived in [25] . It can be interpreted 
as applying adiabatic theory to stochastic systems described by Fokker-Planck equations. Let 
us consider a complex system with hybrid dynamics whose state is fully determined by 

1. N real (continuous) variables x = (xi, ...,#/v) £ ^ N \ 

2. M (discrete) states of a Markov chain. 
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At each time t the state of the system is given by (x, a) G M N x E, where E = {1, M }. The 
dynamics of the MC is given by the specification of the rates, i.e. by specifying the infinitesimal 
generator JC. Therefore the probability of being in state i at time t is determined by the Cauchy 
problem 

f Pitt) = Y^f=iPj(t)K<ji(*)i (7) 

where in general the infinitesimal generator depends on x. In the following we make the basic 
assumption that the rates /Qj(-) are continuously differentiate with respect to each Xi, i = 
1, . . . , N. If we fix a certain state a of the MC, the time evolution of x is determined by a set of 
vector fields f( a \x.) depending on the MC state. Therefore x(t) = (xi(t), ... : x^-(t)) is given by 
solving: 

( d M x-(t) ( \ 

= fi (*i (t), x N (t)), where a = 1, M and z = 1, TV, 

< (8) 

k Xi(0) = x ij0 . 

Here the superscript (-)^ denotes that the time derivative is determined by a vector field de- 
pending on the state a of the MC. We are now in a state to associate an interaction graph to 
the combined hybrid model consisting of both continuous deterministic and discrete stochastic 
variables: 

Definition 3.1. We define the Combined Interaction Graph Xq = (V,E) by 

V{X C ) = V{X D ) U V(X M c) = {xo,X!, . . . ,x N } U {<ti, . . . ,<r M }, 

E(X C ) = E(X D ) U E(X MC ) U E(X D X MC ) U E(X MC Xd), 
ty^A (xi, cr^) G E(X D X MC ) iff ^-JCjj(x.) ^ 0, with i = 0, . . . , N, j = 1, . . . , M, and (a^Xj) e 
E(X MC ^X D ) iff£jfj ai \x)^0, wiihi = 0,...,M andj = 0,...,N. 

This means the combined interaction graph vertex set has as expected the union of the vertex 
sets of the deterministic interaction graph X^ and the Markov chain interaction graph Xmc- The 
edge set of the new hybrid graph Xc also contains the union of the edge sets of Xr> and Xmc, 
but E(Xq) is possibly larger. The set E(X]j — >• Xmc) describes how the continuous variables 
influence the switching rates of the Markov chain, whereas E(Xmc —> %d) describes how states 
of the Markov chain influence the continuous variables. In our examples of Xq (Fig. |ij Fig. [FJ 
Fig. [i| edges contained in E(Xmc —> %d) are drawn by simple broken line arrows, whereas edges 
in E(Xjj — » Xmc) are drawn by bold solid arrows. We also use the convention that if a node in 
V(Xd) influences all the nodes of the MC, then the bold arrow is drawn only once to the box 
surrounding the Markov chain graph. With these concepts we can now conveniently introduce 
the adiabatic or averaging principle: 



Averaging Principle and Averaged Interaction Graph 

Upon the assumption that the evolution of the MC is faster than the deterministic dynamics, 
the macroscopic evolution of the combined system is described by the average dynamics, which 
is given by 
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^ = f> ff (x(i))jf>(x(i)), 
K Xi(0) = xi$. 

Here \i is an invariant measure of /C(x), and i = 1, TV. The invariant measure \i can be a 
convex combination of invariant measures compatible with the initial values {pi(0)}f£ 11 namely 
the support of /i in E is equal to the support of p(0). Note that the average vector field can also 
be written in matrix compact form 

*(f)=/x(x(*))F(x(*)), (10) 
where /i(x) = Oi(x), /i M (x)) and F(x) = {/^(x)}^ 1 ;;;^- 

Definition 3.2. We define the Averaged Interaction Graph X = (V,E) by V(X) = V(Xd) and 
E{X) = E{X D ) U E h with ( Xi , x 3 ) e E u i, j = 1, . . . , N iff 

rl M 

* <T=1 

This essentially means that all (microscopic) interactions of the MC visualised by nodes and 
edges are lumped into edges only by applying the averaging principle, the edges of E^. This is 
precisely illustrated in Fig. [I] The edges in E± are drawn by convention as red (grey) broken line 
arrows, see the examples Fig. [4j Fig. [6] and Fig. [9j 

Remark 3.1. This is an important remark about feedback loops and qualitative behaviour of the 
average dynamics. Because X is again the interaction graph of a deterministic vector field, all the 
known results about the connection between directed loops or cycles contained in X and possible 
restrictions on the attractor sets apply. For example the classical result by Thomas-Soule holds, 
see \TEj . To check the respective conditions it is necessary to look at sign weighted interaction 
graphs, i.e. the existence conditions for arrows need to be refined in the sense that the sign of 
the respective derivatives needs to be evaluated. This is simple in all cases of interaction graphs 
introduced in this paper. Excitation or stimulation symbols in application diagrams correspond 
to arrows with a positive weight, whereas arrows with negative weights correspond to inhibition 
symbols (compare with Fig. [?[). It is presently unclear in which cases the averaging produces 
monotone d ynami cal systems f2$j which generically do not exhibit oscillatory behaviour ( compare 

with section 4-3.1), or when reversely positive feedbacks (positive directed cycles inX) are created 
that are a condition for oscillatory behaviour by the Thomas-Soule theorem. Similar remarks 
could be made for all function classes, for example about the quasi-monotone functions, or the 
ocurence of vector fields allowing chaotic dynamics. In general the sign structure of the Jacobian 
or equivalently the weights ofX are time- and therefore state dependent, with the monotone and 
quasi-monotone vector fields being exceptions. 



3.1 Markov chain structure and the average dynamics 

The MC structure is determined by E and by its infinitesimal generator /C. The state space E 
consists of all possible discrete states corresponding to those molecular (or other) states that 
cannot be described by continuous variables (i.e. concentrations). In the systems of interest 
there are as discussed usually in addition many molecular (or more generally 'discrete state') 
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machines that must be modeled by a discrete space structure independently. Let us now assume 
that there are more than one, say M many such entities. Then each of them will have its own 
state space E^ with i = 1, . . . , M. The total space is therefore the Cartesian product 

E=x i =iE i . (11) 

Each a G E has the form a = (Ji x ...i M = {Oiii •> ■■■> Omz m )i with Oji G Ej. As explained in 
[26] . the infinitesimal generator /C is constructed by assuming that a transition 

a -» a 

involves only one factor space Ej, namely there exits only one index i such that 

G = (Oliu ^i/c-, Omi m ) Cr' = (Oi^, OifeJ, Omi m )- 

This can be called a principle of local effect and is very reasonable for molecular systems. 
Let us consider the simplest case which is E = Ei x E2, where Ej = {0/i,0^} 7 i = 1,2. Set 
^1*2 = (^1*2)? where in order to simplify notation we identify the index with the state. The 
average dynamics then has the following form: 

*(*)= E M( lll2 )(x(t))/ ( ' ll2) (x(t))= £ A*( Ma )(x(t))/ (ili2) (x(t)), 

(iii 2 )GEiXE 2 iiGEi,i 2 GE 2 

where 

^ /i (M2) (x) = l VxGR^. (12) 

(ni 2 )GSixS 2 

Now suppose that the vector field is a linear combination of the form 

f^\ X )=a l9 ^\ X ) + a 2 g^\x) Vx e R N , V(n, i 2 ) £ Si x S 2 , 
with ai, a2 £ R- Then the average vector field turns out to be 



x(t) = aiATi(x(t)) + a 2 ^ 2 (x(t)), 
where the vector fields A\ and A2 are given by 



^( x )= E ( E ^(M 2 )(x)) 5 (ll) (x), 
^(x)= E ( E ^ 2 )(x)) 5 (i2) (x). 



Note that the normalisation ( 12 ) implies that the two expressions 



Mi,(m)(x) = Yl ^1*2) ( x )> 

M2,(i 2 )(x) = ^(M 2 ) ( X ) 



(13) 



(14) 



(15) 



form two probability measures respectively over Ei and E 2 . As discussed in [26] the transpose 
of the generator JC has components that can be written as 
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h,(u>) if 3 =?, 

k {ijWj') = \ k 2,ti3') lf i = ( 16 ) 

if i ^ i' and j 7^ j', 

where fc 1? (.) and fc 2 ,(.) are the entries of the matrices /Ci,/C2 associated to the MC on Si and 
S 2 respectively. 



Proposition 3.1. JJtfte infinitesimal generator JC satisfies relation (16) then the measures /ii(x) 



and /i2(x) are invariant measures of the MC's (/Ci,£i) and (^2,^2), respectively. 
Proof. It is sufficient to prove the result for one MC. Let us compute 

This is equal to 

Yl fe (i,iiii)M(ii,i 2 )- 
The last expression can be rewritten as 

E E E hiMi'^i^A)- (17) 

^2^X12 *2^ S 2 

Now using the definition (TTgI) , expression (I17l becomes 



J2 hni[):(i2i f 2 )^i[A) = °- ( 18 ) 

iiGEi i 2 GS 2 z' 2 e£ 2 

The latter equation is satisfied because /i is an invariant measure on (/C, Si x E2), and therefore 

□ 

The result of this paragraph can be easily generalised to systems where the MC is constructed 
by the product of M MC's. We leave this extension to the reader. 

4 Examples 

This section is devoted to different examples of increasing difficulty which will explain and demon- 
strate the approach and introduce examples of interaction graphs. The first two examples con- 
tains the analysis of a so-called feed-forward-loop (FFL) whose modelling has been introduced 
in [20]. In this context we look at both an AND and OR gate. We took these two examples to 
show that the average dynamics is a relatively natural way to attach testable realistic dynamics 
to a well known genetic wiring diagram, or mathematically speaking, to a directed graph. Our 
final example will present a new point of view on modelling a genetic oscillator. In our models 
each single gene is described through a set of finite discrete states with the dynamics given by 
a Markov chain. The switching of genes is triggered by molecules described as a continuum, 



15 



typically interpreted as transcription factors that are sufficiently abundant to be described by 
a continuum variable. In this setting we show that the average dynamics for the clock allows 
to demonstrate circadian rhythms using a two-gene system - without the need to introduce the 
assumption that genes are modelled as 'gene concentrations'. This is unfortunately an underlying 
assumption for most existing models found in the literature. Note that in the modeling step we 
usually interpret transitions as reaction schemes, involving classical mass-action reaction for con- 
tinuous variables, stochastic transitions for MCs, but also combined reactions incorporating both 
mass-action kinetics and stochastic transitions. See [18], [25] for a more detailed introduction to 
this extended reaction scheme formalism. 



4.1 Feed- forward loops 

In [20 , Mangan and Alon analysed so called feed-forward loops (FFL) describing genetic circuits. 
We interpret them here as a single wiring diagram only, assumed to be independent from any 
other process occurring in a cell. This means we are not interpreting them as a small subgraphs 
of a larger graph, i.e. as so-called motifs. Such loops if isolated in this sense are generically 
described through diagrams like Figure [2] In principle the choice of the associated dynamics is 
open, see the discussion in [18 . One can first select the interaction graph, then the dynamics, 
or vice versa. There is in other words no one-to-one relationship between the graphical model 
and the equation determining the time evolution. A better interpretation of the interaction 
graph is that the feed-forward loop is determining a possible class of equations. The question 
is then whether this class has distinctive common features which can be determined - or not. 
The generic class of equations associated to the FFL in the literature is chosen as a system 
of ordinary differential equations which associated vector field is compatible with the graphical 
representation of the FFL. In [20 the diagram in Figgis associated to a dynamics given by 

y = by+ Py /(X*, k X y) Oiy i/ , 

i = b z + j3 z GO*, k xz , y*, k yz ) - a z z, ^ J,) 

v 

where 

x if Sx = 1, 
if Sx = 0. 



The 2/*, the transcription factor Sy and b yi b Zl f3 y , f3 z ,a y and a z are all positive constants. 

i-linear 

f(u,k) 



The key elements of (19) are the non-linear functions / and G. For an activator we assume 

(u/k) H 



(l + (u/k) H[ 
whereas for a repressor the non-linearity becomes 

f{U ' k)= (l + (u/k)f 

The gate function for an AND-gate is given by 



G = f(x*,k xz )f(x*,k yz ). (20) 

If there are two transcription factors competing for binding to the promoter region, a possible 
choice for an activator is 
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OR gate 





Figure 2: AND gate and OR gate according to U. Alon et al. 



and for a repressor 



(u/k u ) H 



l + (u/k) H + (v/k) H ' 



f c {u] k Uj k Vj v) 



l + (u/k) H + (v/k) H ' 
For the OR-gate in [20] the G transfer function is 

G = f c (x* y*) + fc(x*-,k 

yz i &xz i y*)- ( 21 ) 
Next we demonstrate that the averaging principle can be used to derive effective dynamics with 



qualitative properties similar to (20) and (21) giving a new interpretation and justification to 



the choice of non-linearities selected to describe the dynamics. The derivation is based on on a 
sequence of assumptions necessary to interpret graphs like Fig. [2] The interpretation is essentially 
related to the description of the processes at a semi-microscopic level where gene activation is 
modeled. Because assumptions are made on the microscopic scale, they are open to experimental 
(molecular) validation. 



OR-gate 

We assume that X and Y are competing for a single binding site P. This scenario can be 
modelled by the following set of reactions: 



X + P ^ PX PX + z 



hi 



r + P ^ PY PY + Z 

h 2 

Z 
-+ h * Z 



(22) 



The OR-gate reactions can be described using the graph in Fig. [3j Note that in this case we 
need to include an environment vertex. In the graphical representation the environment vertex 
includes also the species X and Y which affect the MC but do not have an explicit dynamics. 
The MC has the state space 



X = {P,PX,PY}, 
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and the infinitesimal generator is given by 



-kix - k 2 y hx k 2 y 
JC= | h x -h x 

h 2 -h 2 




Figure 3: The OR gate combined interaction graph Xq- Note that the environment node 
(x, y) has two components. 

The deterministic dynamics is chosen as simple as possible: 

z(t) = b z — S z if the MC is in state P, 

z(t) = b z - 5 z + p x if the MC is in state PX, 

z(t) = b z - S z + P y if the MC is in state PY. 

The invariant measure of the MC is 



H(x,y) 



hxh 2 



kxh 2 x 



k 2 hx y 



hxh 2 + k x h 2 x + k 2 hxy' h x h 2 + feift 2 a: + k 2 h x y' h x h 2 + feift 2 + &2^i 2/. 
and the average dynamics turns out to be 



AiiAi2 + kxxh 2 + «22//ii 
The transfer function of the OR-gate then becomes 

kxxh 2 



Py 



G = P X 



hih 2 + k\xh 2 + k 2 yhi 



k 2 yhx 

hxh 2 + kxxh 2 + A:22/^i 

k 2 yhx 



hxh 2 + &ix/i 2 + fe 2 2/fei 



The average dynamics has an associated interaction graph X as visualised in Fig. [4] 
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Average deterministic dynamics 



Figure 4: In the interaction graph X associated to the average dynamics the whole MC is 
substituted by a second (henceforth new) edge directed the environment state to Z. 

AND-gate 

In this case we assume that both X and Y can temporarily bind to site P. This can be modelled 
with the following set of reactions: 



X regulates Y: 
X 



ki 

P^ PX 



PX -> 6 « PX + Y 



X and Y regulate Z: 



Y 
X 
Y 



PY 



fci 



- PY ^ PXY 

hi 

PX ^ PXY 

h 2 

PXY PXY + Z 



(23) 



Degradation of Y and Z: 
Z -^ 5z 
Y -^ s y 



Basal expression of Y and Z: 
-+ bz Z 
-^ b y Y 



(24) 



These reactions can be again described through the interaction graph given in Figure m Also m 
this case an environment vertex has to be introduced. The Markov chain has state space 

E = {P,PA, PY.PXY}. 
The corresponding infinitesimal generator is given by 



AC: 



V 



k\x — k2y 


k\x 


k 2 y 







-hi - k 2 y 





k 2 y 


h 2 





—h 2 — k\x 


k\X 





h 2 


hi 


-hi - h 2 



The deterministic dynamics is again chosen as simple as possible: 
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y(t) = b y - S y y if the MC is the state P, PY and PXY, 
y(t) = by - 5 y y + /3 y if the MC is the state PX, 

= 6* - 5 Z z if the MC is the state P, PX and PY, 
z(t) = b z - 5 z z + j3 z \i the MC is the state PXY. 




Figure 5: The AND gate combined interaction graph Xq. 
The invariant measure of the MC is 



fi(x,y) 



h 2 hi 



k\xh 2 



h 2 h 1 + k 1 xh 2 + k 2 yh 1 + k x xk 2 y" h 2 hi + fcix/i 2 + &22/fti + hxk 2 y' 
k 2 yhx k\xk 2 y 



h 2 h\ + k\xh 2 + k 2 yh\ + k\xk 2 y^ h 2 h\ + k\xh 2 + k 2 yh\ + k\xk 2 y ) ' 
and the average dynamics turns out to be 

k\xh 2 



y = b y - 5 y y + f3 y 
z = b z - 5 Z z + f3 2 



h 2 h\ + k\xh 2 + ^22/^1 + kixk 2 y' 
k\xk 2 y 



h 2 hx + &i£/i 2 + ^22/^1 + hxk 2 y 
This implies we obtain an AND-gate with transfer function 

G = ^ fcia;fc22/ 

/i2^i + k\xh 2 + k 2 yh\ + k\xk 2 y 

The averaged dynamics has an associated interaction graph X as given in Fig. |g] 
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Figure 6: The averaged interaction graph X associated to the average dynamics of the AND 
gate. Two new edges are formed between Y and Z. These links correspond to new effective 
interactions as described by the edge set E±. 

4.2 An application to genetic oscillators 

Cyclic behaviours are important elements for the understanding of regulation inside any biological 
process. It is no surprise that therefore oscillations are at the centre of many different approaches 
to their modelling and analysis (see [2j El [7] [12] [TJ [30] [23] E]). We draw particular attention 
to [30] where a genetic oscillator using a system of two types of interacting genes has been 
designed. The structure of the model is presented in Fig. [7] The two genes gene a and geneR 
express two proteins A (activator) and R (repressor). The activator A enhances gene a and 
geneR. The activator A inhibits the repressor R by forming the complex C. In [30] the model is 
constructed by assuming that there exist two uniform populations of genes and thus mass action 
kinetics is used to describe the activation process. We now show that the unnecessary underlying 
assumption of a uniform continuous population of genes in the cell can be overcome by using 
the average dynamics approach. In fact we re-formulate the model by considering just two genes 
described by a MC on a finite state space. 




Figure 7: Genetic Oscillator Model by N. Barkai et al. 
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4.2.1 The associated reaction scheme 

We keep the same notation as in [30] in order to simplify the comparison. So we use the following 
notation: 

• A is the activator protein and Ma its corresponding mRNA, 

• R is the repressor protein and Mr its corresponding mRNA, 

• C is a complex formed by A and R. 

Each gene can be either active or inactive, Da, Dr denote the inactive states, and D' A , D' R the 
active states. We now collect all necessary reactions to model this situation: 

Gene activation: 

A + D A ^D' A , 
A + D R ^D' R 

Transcription: 

D A ^ aA Da + Ma, 
D r ^ r D R + Mr, 

Translation: 

M A ^ a A, 

Regulation and inhibition of protein A: 

A + R -> lc C, 

Degradation: 

M A ~^ 5ma 0, 
A -^ Sa 0, R -^ Sr 0. 

We next need to introduce the following concentrations for smaller molecules in the system 
and their complexes: 

a=[A], c=[C], r = [R], m A = [M A ], m R = [M R }. 

4.2.2 A combined finite and infinite state space formulation 

To apply the averaging principle it is necessary to identify finite and infinite variables or degrees of 
freedom. As discussed for genes interpreted as large molecular machines in perhaps single copy 
number the mass action kinetics typically cannot be used to describe the activation process. 
Therefore it is a natural choice that the two genes gene a and geneR form the finite degrees 
of freedom whose time evolution will be given by a finite Markov chain. The concentrations 
of A,C,R and Ma, Mr will be naturally the continuous variables as these molecules are much 
smaller and very abundant relative to the genes. We first describe the MC structure, and then 



D' A ^ D' A + M A 
D f R ^D> R + M R 

Mr -+P r R 

c -> 5a r 

Mr ~^ 5mr 0, 
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state the equations for the concentrations a, c, r derived by applying the principles of mass action 
kinetics. The structure of the dynamics can be represented through the combined interaction 
graph in Figure [8] The structure of the links will be made clear by the construction of the MC 
and of the related deterministic dynamics. 



/1m!drV- ♦fD'A.DR 



ii \ / Ji 









MA I 
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D 


I R J 

( C \ 
sterrr in st c cyanics 



Figure 8: Combined interaction graph Xc for the VKBL model. 



4.2.3 MC structure 

The Markov chain of the VKBL model has a state space given by all possible states of the couple 
(gene a , gene r). gene a and genen have state spaces respectively 

Z A = {D A ,D' A }, E R = {D R ,D' R }. 

Therefore the total space is equal to 

S = SixS fi = {{D A , D r ), (D' a , D r ), (D a , D' r ), (D' a , D' r )}. 

The reactions associated to gene activation determine the infinitesimal rates of the MC. We 
consider transitions in which either Da or Dr is affected as illustrated in Fig. [8] Transitions 
like (Da,Dr) — >• (D' A ,D' R ) are excluded because they occur with a very small probability. In 
fact it is very unlikely that two molecules A bind simultaneously on the two different existing 
binding sites. In order to appreciate the effects of describing A protein with its concentration a 
let us take the matrix K7 (the transpose of the infinitesimal generator) written in terms of the 
A copy number always denoted by a: 



JC 1 



( -«(7i? + 7A) 




A 

-6 A - (a - 1) 7,r 

- l)7i* 






Or 
-9a - Or 
Oa 



Or \ 


- 1) 1A 

9 R - (a - l)j A J 



(25) 
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Remark 4.1. The copy number a is an integer greater or equal 1. Note that we are not applying 
the average principle here in order to be able to discuss the transition of small particle to large 
particle copy numbers. 

The invariant measure is given by 

, x _ f OaOr ORjAa 

\OaOr + 6 R ^ A a + a 2 7^7A - jR^Aa + Gacl^r ' 6 A 6 R + R ^ A a + a 2 ^ R ^ A - 7h7a« + Oaohr ' 

ajR (a - 1)ta 0AQ7H ^ 

OaOr + Qr^aol + a 2 7 jR 7A - iR^Aa + Gacl^r ' #a#,r + R ^ A a + a 2 j R j A - 7#7a« + Qacl^r J 

(26) 



4.2 .4 Large values of A concentration and independent genes 



Let us now consider the regime of large copy number of A, namely a >> 1, this implies that A is 
well described by concentration. In this case we could approximate a — 1 c± a. This would yield 
a MC with a new /C T given by: 



JC 



fT 



( -a(7i? + 7A) 

J Ad 





V 



a7i? 



pa 

-6 A - a<y R 






Or 
-@A — Or 

Oa 



vr 




\ 



CIJA 

-0 R - aj A ) 



(27) 



JC /T describes a product of two independent chains. In fact let us recall the structure of the 
spaces Yia = {Da,D' a }, = {Dr,D' r }. The matrix JC fT can be rewritten in the form (16) 
using the the following matrices 



K7 A 




Oa 



and 



In the regime a 



fi(a) 



I -aj R 6r \ 

V a lR -Or J 
a the invariant measure is equal to 

Oa Or a >y A Or 



(28) 



(29) 



+ o7a)(0j* + uIrY {0 a + cij A )(0 R + a^ R )' 



or 7 'a 1r 



aiR Oa 



li A {a) 



aj A )(0 R + a j R ) ' (0 A + a ja)(0r + a j R ) 

T 
B 

CI J A 



Now note that matrices IC A and /C^ have invariant measures which are given by 



A 



Oa + clja Oa + clja 



and /jLr(o) 



cilR 



0R + a<y R R + a<y R 



One can easily verify that the components of the invariant measure fi(a) satisfy (15) and factor 
into the components of iia{o) and fiR(a): 
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Ma) = (ti\a)^(a) , f,f(a)^(a) , $\a)$\a) , ^\a)^\a)). (30) 

This corresponds to have a gene activation formed by the product of two independent MC's 
whose generators are respectively ICa and IC R . Intuitively we can say that in the regime a » 1 
the molecules activating the two genes are so abundant that genes do not compete and behave 
independently. The mechanism can be represented by the two reactions 

A + Da D' A , and A + D R ^ D' R . (31) 

0a Or 

Remark 4.2. Note that if a » 1 is not valid, then the dependence of the MC's can be interpreted 
as an cooperative effect, in the sense that the two MC's interact through the presence of the A 
molecules. 



The product of probabilities in (30) indicates that the two MC's are independent whenever 



the number of A molecules is so large that the reactions (31) can be considered statistically 
independent. In general one can envisage the possibility of having systems with large MC s 
interacting with a certain type of molecule. The chains do obviously compete for such molecules. 
As long as the number of molecules is small then the chains are in competition. But if the number 
of particles becomes large then the chains are effectively de-coupled. The stochastic analysis and 
the related algebra of this problem will be developed in a different paper. 



4.2.5 Deterministic dynamics 

In the regime a » 1 the dynamics of A, R and C is obtained by using mass action kinetics. It 
turns out that the dynamical equations are affected by the MC only through the states Ma and 
Mr (modelling the transcription process), and they are given by 

a(t) = -S A a(t) - 7c a(t) r(t) + Pa m A (t), 



c(t) = -6 A c(t)+7 C a(t)r(t), (32) 

r(t) = -S R r(t) - 7c a(t) r(t) + f3 R m R (t) + S A c(t). 

Note that since the equation for a(t),c(t) and r(t) do not explicitly depend on the MC state, they 
are not affected by the averaging procedure. Of course the dynamics of the mRNA is dependent 
on the state of the genes. Let f^ be the vector field corresponding to the state a G S. The 
transcription processes are governed by 



m A (i) 




CtA 


- S M a m A (t) 


- p A m A (i), 


rh A (t) 


= f(D> A ,-) 


= a' A 


- S M a m A (t) 


- p A m A (t), 


m R (t) 


= f(,D R) 


= OLR 


- 5 M r m R (t) 


- Pr m R (t), 


m R (t) 




= a R 


- $mr m R {t) 


- Pr m R (t). 



4.3 Average dynamics for the VKBL model 

The average procedure does affect only tua and m R . Indeed we find that 
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*M*) = V(D A ,D)f 
De{D R ,D' R } 



(D A ,D) 



(D f A ,D) 



and 



m R (t) 



Er(D,D R ) i f| r(D,D') 

De{D A ,D' A } 

After some algebra the average dynamics can be shown to be 



a(t) 

C(t): 

m 

m A - 
rh R - 



-S A a(t) - 7c a(i) r(t) + fi A m A (t), 
-5 A c(t)+'yca(t)r(t), 

Sr r(t) - 7c a(t) r(t) + p R m R (t) + S A c(t), 
a A A a f A j A a (t) 



Oa + 7a a 0) 0a + 7a a (*) 
olrOr a' R a {t) j R 



(5ma + fi A )m A (t) , 
(^Mi? + Pr) m R (t) . 



(34) 



Or + a (t) 7^ + a (t) 7^ 

The latter set of equations corresponds to the equations as stated in [30], giving them a new 
interpretation. 
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Averaqe Dynamics 



Figure 9: The averaged interaction graph X for the VKBL model. 



By inspection of (|34j) the averaged interaction graph X can be constructed and becomes as 
shown in Fig. [9j One can verify that new links are present due to the effective transcription 
rates depending on a. 

4.3.1 Numerical analysis of the model 

The genetic oscillator described by the average dynamics can be investigated numerically. We 
use the values of parameters given in [30] and make the following choices: 
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Figure 10: A typical average dynamics time evolution for proteins A (red line 1), C (blue line 2) 
and R (green line 3) following a limit cycle. 



We integrate the average dynamics (34) with an ODE solver numerically (using GNU Octave). 
Here we have chosen the initial conditions 



a(0) = 0, c(0)=0, r(0) = 0, m A (0) = 0, m R (0) = 0. 

In a second step we used Dizzy to simulate the average dynamics using Gillespie's algorithm, see 
Fig. [TT] This means we interpret the average dynamics vector field as chemical reactions where 
the particle number involved in the reaction stays finite. Similarly to the model presented in [30] 



we also found in this case numerically a limit cycle. This is shown in Fig. 12 In fact even though 
protein abundances now fluctuate in time this limit cycle is very close in average norm to the one 
created by the average dynamics. It holds of course that this similarity is higher if more particles 
are used in Gillespie's algorithm to simulate the dynamics. Note also that this limit cycle can 
be obtained analytically from the average dynamics directly when using the Hopf bifurcation 
theorem. One can verify that for the given parameter values there is one single unstable steady 
state which is located in the interior of the 'stochastic' limit cycle projected into the c, r-plane. 



5 Conclusions 

In this paper we presented the average dynamics as a framework to derive macroscopic models 
based on typical microscopic assumptions often occurring in complex systems theory. The system 
is separated into finite-state machines (themselves finite in number, typically occurring in a single 
copy number only), and system components that are associated with very abundant particles 
moving randomly and spatially homogeneously. These latter particles are given as concentrations, 
and they are typically describing particles communicating between the finite-state machines, 
creating switches between the states. The average dynamics is then naturally constructed out 
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Figure 11: A typical random time evolution of proteins C and R using Gillespie's algorithm in 
a stochastic version of the average dynamics (Randomness is due to finite particle numbers, not 
because of stochastic transitions in the original Markov chain part of the model). 




Figure 12: Projection of the 'stochastic' limit cycle into the c, r plane. 
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of a collection of deterministic vector fields and depending on a Markov chain steady state. 
The construction is a very useful tool to derive effective models that approximate stochastic 
systems described by particular reaction schemes acting on a mesoscopic or microscopic scale. 
One of the best examples of such a situation is the use of combinations of Markov chains to 
describe conformational changes of large bio- molecules. In this case the average dynamics allows 
to establish a very modular approach to study complex bio-molecular diagrams. Such models can 
be tested with data derived from measurements on different scales, and are therefore much better 
to test empirically. The framework described can also be interpreted as allowing to assign more 
or less automatically deterministic dynamics to given interaction diagrams, where the modeler 
can make choices whether a variable can have a continuous range, or a discrete range, [T5l [T6]. 
The diagram then has a mathematical precise interpretation as an interaction graph associated 
to the dynamics. It can be further studied to investigate the qualitative behaviour of the system. 
The interaction graph is a more precise tool to identify feedback loops when compared to the 
motifs which as a static (i.e. time-independent) concept are sometimes inappropriately used in 
the same or similar role in Systems Biology [18 j. We made the step to use such graphs also in a 
multi-scale and hybrid dynamics setting. 
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